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Abstract: We develop a variational wave function for the ground 
state of a one-dimensional bosonic lattice gas. The variational theory is 
initally developed for the quantum rotor model and later on extended 
to the Bosc-Hubbard model. This theory is compared with quasi-exact 
numerical results obtained by Density Matrix Rcnormalization Group 
(DMRG) studies and with results from other analytical approximations. 
Our approach accurately gives local properties for strong and weak 
interactions, and it also describes the crossover from the superfluid 
phase to the Mott-insulator phase. 
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1 Introduction 

During the last years a spectacular development in the storage and manipulation of 
cold atoms in optical lattices [1, 2] has taken place. Greiner et al. [1], to name one 
important example, succeeded in experimentally driving a quantum phase transition 
between a superfluid and a Mott-insulating phase in bosonic systems. This experimen- 
tal progress has revived the interest in the Bose-Hubbard model [Eq. (1)] as a generic 
Hamiltonian for strongly correlated bosons, by which the quantum phase transition can 
be described [3] . The Bose-Hubbard Hamiltonian has been used previously in condensed 
matter physics to study the adsorption of noble gases in nanotubes [4] , or Cooper pairs 
in superconducting films with strong charging effects [5, 6]. In this context a lot of work 
has already been done to characterize the quantum phase transition, the statistics, and 
the low-energy excitations of the Bose-Hubbbard model [7, 8]. However, new interesting 
questions arise now due to the good tunability of the experiments with optical lattices. 
In particular, it becomes possible to study time-dependent processes such as driven 
quantum phase transitions [1]. A theoretical understanding of such phenomena is chal- 
lenging, since the characteristics of the superfluid phase — where the atoms tend to 
delocalize throughout the lattice and large fluctuations in the local density exist — , 
and the Mott-insulating phase — where the number fluctuations decrease, and a gap 
in the excitation spectrum opens - must be covered at the same time. Both regions are 
separated by a non-analyticity of the spectrum, which implies that a perturbative study 
[9, 10] works best in strong coupling limit, while a Hartree-Fock-Bogoliubov mean field 
works best in the superfluid regime. In addition it is possible to develop a mean field 
theory [11, 12] based on a Gutzwiller ansatz [3]: this reproduces the mean field theory 
in the superfluid limit, as well as the limit of infinite interaction, which raises the hope 



that the theory also interpolates properly between these limits. 

In this paper we develop a variational description of the ground state of an ensemble 
of cold atoms in an one-dimensional optical lattice. Our trial wavefunction treats the 
connections between neighboring sites as entities which decouple in the limit of infinitely 
large lattices. We apply this technique first to the quantum rotor model — which de- 
scribes the lattice for large and commensurate occupation per site [5] , and has also been 
used to describe an array of Josephson junctions [13] — , and next to the Bosc-Hubbard 
Hamiltonian. The accuracy of the variational theory in both models is confirmed by 
several comparisons. In the quantum rotor case, we use a spin wave approximation in 
the limit of weak interaction, and a first order perturbation theory in the limit of weak 
tunneling. For the Bose-Hubbard model we compare against results obtained apply- 
ing the quasi-exact, numerical DMRG method to one-dimensional lattices with up to 
128 sites, and also with calculations based on the Gutzwillcr ansatz. Our conclusion is 
that the variational picture of self-regulated connections between sites provides a rather 
cheap and simple, but very good description of the local properties of the system in the 
supcrfluid and insulator regimes, and a fairly good interpolation across the quantum 
phase transition. It cannot describe, however, the algebraic decay with distance of the 
off-diagonal elements of the one-particle density matrix. The simplicity of the method 
suggests a possible generalization to higher dimensionalities and other physical models. 

The outline of the paper is as follows: In Sec. 2 we introduce the quantum rotor 
model as a possible limit of the Bose-Hubbard Hamiltonian. Next, information about 
the ground state of the quantum rotor model is obtained variationally as the solution 
of a Mathicu equation. We can estimate energies, correlation functions and length, and 
the variance of the density as a function of the only free parameter. A comparison with 
perturbative estimates demonstrates the accuracy of the method when computing local 
properties. Since the quantum rotor model is only an approximate description of the 
optical lattice, in Sec. 3 we develop a similar variational theory for the Bose-Hubbard 
Hamiltonian. After bringing the Bosc-Hubbard Hamiltonian to an appropriate form, 
we can estimate the local properties of its ground state. The variational solutions are 
compared in Sec. 3.3 with the results of DMRG studies of the Bose-Hubbard model. 
We confirm that the variational method describes very well the local properties of both 
the Mott-insulator and the superfluid regime, and provides a fairly good interpolation 
across the phase transition. Finally, in Sec. 4 we summarize our results and comment 
on possible extensions. 

2 Quantum phase model 

2.1 Relation to the Bose-Hubbard model 

In this section we show the equivalence of the Bosc-Hubbard model 



and the quantum rotor model for large and integer occupation n [5] . In our notation, M 
is the number of lattice sites and N = nM the number of atoms. Both the Bosc-Hubbard 
model and the quantum rotor model show a phase transition due to the interplay be- 
tween the kinetic term proportional to J and the interaction term proportional to U . 
For convenience we have subtracted the ground state energy in the perfect insulator 
limit U / J — > oo. 
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If we expand a configuration of the lattice using Fock states 
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and the number of particles per lattice site is large, > 1, we may approximate the 
hopping terms as follows 
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Here, A^n) = \n ± 1) are ladder operators and V projects on states with non-negative 
occupation numbers, > 0. To lowest order the error |A;j) is 
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and its norm is bound by 
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where oi = ((n; — n) 2 ) is the variance in the number of particles per lattice site. For 
the approximation (3) to be valid, the uncertainty in the number of atoms must be 
small compared to the mean value, n > <r, and the interaction energy must exceed the 
neglected terms, XJn{n — 1) » J a. 

Following the previous procedure the Bose-Hubbard model becomes 



H QR = v E -2J P (A+ +1 AJ + A+Aj +1 ) 



(6) 



Here p = \/n(n + 1) is approximately the density, A| = OjCij — n is essentially the 
number operator and we have used that ■ A^\%\i) =0 when we work with states that 
have a fixed, commensurate number of particles. In the following we will define the 
energy per lattice site as 

e=±{H QR ). (7) 

Since the physically interesting states will be concentrated around large occupations, 
nk = n, the usual step now is to drop the projector, V, and move to the basis of phase 
states, defined by 
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In doing so, we obtain the identification Ad- 
duces the usual representation of the quantum rotor model 



H, 



QR 



E 



-2Jp cos(cj)j 



u d 2 

2 d(j> 2 



(8) 

-id/d(f>j, which pro- 



(9) 



with the associated state writen as 

\i>) = (2tt)- m / 2 J d M 4>e in E0fc *(<?)|<?). 
A similar derivation is possible using path integrals [15]. 
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Fig. 1. Instead of working directly with the population of each well, n k , we can use 
other quantum numbers, w^, defined by the relation = — «»fe_i + n, and 
which behave like a set of chemical potentials acting on the barriers that connect 
neigboring sites. 

Variational ansatz 



In this section we estimate the properties of the ground state of Hqr variationally. 
Due to the previous splitting (10), any wavcfunction ty(<f>) can only depend on the 
phase difference between neighboring wells, £j = 4>j+i — 4>j. These new quantum vari- 
ables describe the connections between neigboring sites. In the limit of large lattices 
it seems reasonable to assume that these connections become independent from each 
other adopting the product state 
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This representation becomes exact in the Mott-insulating limit, U/J — > oo, where 
hmottiO = lj an d in the superfluid limit, U / J — > 0, where h s f{^) = X)nGZ<^(£ — 27m), 
as can be verified by direct substitution in Eq. (10). 

Even though the phase representation is the best one to find a trial wavefunction, 
it is not the optimal one for performing computations. It is instead more convenient to 
work with the variables which are conjugate to the phase differences £j = 
These are the new quantum numbers, Wk, given by 



Pj+i 



n k = w k -Wk-i+n. (12) 
In terms of these numbers, the ansatz (11) reveals itself as a simple product wavefunction 
\i>) = \h)^ M ^ =Y / h Wl --- h WM _, K) ® • • ■ ® \wm-i), (13) 



with coefficients given by the Fourier transform 
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As sketched in Fig. 1, the Wk play the roles of chemical potentials which are established 
between different wells: the difference between the potentials on the extremes of a site 
gives the fluctuations over the mean and commesurate occupation h. In this picture 
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where £- 



\w ± 1) are new infinite-dimensional ladder operators and E z |w) 



By minimizing the energy associated with Hqr over all states within a given ansatz 
we can both obtain an upper bound to the energy of the ground state and approximate 



its wave function. A simple computation with our product ansatz leads to the result 
[compare (7)] 

e[h] ^ -4JpRc(£+) + t/((S 2 ) 2 ) - U{T, Z ) 2 , (16) 

where the expected values are computed over a single connection, (E 2 ) = w\h w \ 2 , 
and the wavefunctions are assumed to be normalized, ^2 W \h w \ 2 = 1. Since the stationary 
states have a well defined parity, ht- w \ = (— l) p h w , the optimal variational state must 
satisfy the linear equation 

-2Jp(h j+1 + hj-i) + Uj% = Eesthj, (17) 

which is nothing but the Fourier transform of a Mathieu equation 

U d 2 



2 ae 
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h{i) = e est h(£), (18) 



The estimated ground state energy per site is given by the lowest eigenvalue of cither 
equation. 

Using the product ansatz and the same approximations required to derive Hqr, we 
can also compute other properties of the ground state. For instance, the variance of the 
number of atoms per lattice site 

a 2 = ((a]a,-n) 2 )=2((£ 2 ) 2 ). (19) 

and the correlation functions 

(aj+i%-> = p(Vj)=m, (20) 

1 3+1 

(21) 




which decay exponentially with the distance. This implies that the ansatz (13) only 
describes properly the decay of the correlations in the Mott-insulating regime, since the 
correlations in the superfluid regime follow a power law decay. However, as we will see 
below, local magnitudes (a, jx, e . . .) are properly estimated even if long-range ones are 
not. 

We have solved Eq. (17) numerically in a truncated space. The results are summa- 
rized in Fig. 2, where we also plot reference estimates arising from two other analytical 
methods. In the limit {/ > J we compare with a first order perturbative calculation 
around the solution = 1, which is possible thanks to the energy gap of order 0(U) 
in the excitation spectrum. In the limit J > [/ wc rather use a spin wave or har- 
monic approximation in which the cosine-term of the Hamiltonian Hqr is expanded up 
to second order in the phase difference between neighboring sites [see Sec. 2.3]. This 
approximation is valid in the superfluid regime, where the phase does not vary much 
between neighboring wells. From the graphical comparison we see that, as expected, the 
variational wavefunction provides a fairly accurate description of the ground state of 
the quantum rotor model in both the superfluid and insulating limits. As a side note, 
we must remark that this ground state has a divergent fluctuation of the number of par- 
ticles per site as 3 — > 0. This disagrees from the expected behavior of the ground state 
of the original Bose-Hubard Hamiltonian, and it reminds us that Hqr can only model 
the atomic ensemble when the variance, a, is small compared to the mean occupation 
number, n. 
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Fig. 2. Estimates for (a) energy energy per lattice site and (b) density fluctuations 
of the quantum rotor Hamiltonian (6) obtained with the variational method (solid), 
and perturbative calculations for U <g J (dashed) and U 2> J (dots). 

2.3 Harmonic approximations to the quantum phase model 

In the limit J 3> U it is possible to estimate the ground state of the rotor model 
(6) analytically. Since we are deep in the superfluid regime, the wavefunction will be 
concentrated around the line <pi — (f>2 — ■ ■ ■ = (f>M , and we can approximate 
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If we remove the periodic boundary conditions on 
Hamiltonian may be diagonalized, H = ^2 k &k 



and change variables, the preceding 
1 with frequencies given by 
1< 2k < M which 



A\A k 



u>k = y/8JpU | sin(7rfc/M)|, where k is an integer in the range —M 
labels the different values of the momentum in the lattice. The ground state energy [Fig. 
4(a)] may be estimated as the zero-point energy of the harmonic oscillator. For large 
M, the sum over k may be replaced with an integral, giving 
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The variance of the number of particles is related to the expectation value of the mo- 
mentum using the same procedure as above 
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The Bose-Hubbard model 



In this section we apply to the Bosc-Hubbard Hamiltonian the techniques that were de- 
veloped in Sec. 2. We will do it in three steps: First we will develop a phase representation 
which is valid for all occupation numbers. Next we will prove that this representation is 
equivalent to a similarity transformation of the Hamiltonian which brings it to a form 
similar to (6), at the price of losing Hermiticity. Finally we will show how to implement 
the ansatz of independent connections (13) to produce estimates for the usual set of 
observables (e, 71, a), which are to be validated with DMRG calculations. 



3.1 Coherent states 



The phase coherent states \<p) are defined by 

(n\<f>) = e in +ly/n\. (25) 



Unlike the phase states defined in Sec. 2, they are not orthogonal to each other, 
exp [e 1 ^ - ^- 1 ] , but they form a complete basis, so that an expansion like (10) is still possi- 
ble. A nice property of the coherent states is that we can rewrite the operators a, at, a' a, 
etc, in terms of the phases very easily. For instance, a\<j>) = e l< ^\(j)), cfi\<f)) = — ie~ l ^^|</>), 
and a) a\<j)) = —i-J^\(f)). Using this representation, we obtain an effective Hamiltonian for 

the wavefunction i.e. H B hH) = {2ir)- M / 2 j d M 4>e lfl ^ 'M#* oh *(<?)]|<?), which is 

of the form 



Hcoh — J 
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Here H l coh stands for the tranpose of H co h- This operator was already used in Ref. [28] 
to study the Bose-Hubbard model with only two sites. On the one hand, it is a non- 
Hermitian operator 1 and we cannot do a simple variational study. On the other hand 
the Hamiltonian still depends on the phase differences, and it is reasonable to look for 
approximate cigcnstates which have the form (11). This will be done in the following 
section. However, since working with phase variables is inconvenient, we will develop a 
representation similar to that of Eq. (15) in the following. 

3.2 Variational procedure for non-Hermitian operators 

In this section we will find the best variational function which has the product form of 
Eq. (11). However, as it happened in Sec. 2, instead of working with phase variables 
it is more convenient develop a representation of the Bose-Hubbard Hamiltonian in 
terms of connections. This is once more a two-steps process. First we use a similarity 
transformation suggested by the definition of the coherent states 




(27) 

Since Oa.O- 1 = At and OaJO- 1 = A z -Af, we find 

J J J J 3 

Hcoh = OHbhO- 1 (28) 
= -j£(^+n)4 + A- + !]T^) 2 . 

The Hamiltonians (26) and (28) are equivalent: while one is defined in terms of phase 
variables, the other one is defined with occupations numbers, and both are related 
by a Fourier transform. The second and final step is to rewrite everything in terms of 
connections, using the quantum numbers from Eq. (12), and the relations = 
E y = z(I]~ — S + ). The result is a decomposition of the Hamiltonian 

Hcoh = H!+H 2 , (29) 



lr The hermiticity of Hbh is maintained due to an implicit projection that takes place when we 
reconstruct the state H BH \tp) from Hl oh ^($) (See Ref. [28]). 



3 

into terms which are local, Ifi, and terms which involve pairs of connections, Hi- 

For the quantum rotor model we proved that the optimal product wavefunction 
was an eigenstate of a Hamiltonian which did not couple connections, like H±. The 
difference now is that, since the operator H c0 % is not Hermitian, we cannot establish a 
variational principle and that proof is no longer valid. Nevertheless, we will again propose 
a variational ansatz which is an eigenstate of the local operator H\\h)® M = Me est \h)® N . 
Using the following equality 

. (#gggW • (x\0- 2 H coh \x) 
£ = mm ——- = mm — — , (30) 

||?Mj 2 x#o (x|o^ 2 |x) 

and the product ansatz \\) = \h)® M , we arrive to an upper bound for the lowest 
eigenvalue of the Bosc-Hubbard Hamiltonian, expressed in terms of the non-Hcrmitian 
one 

eo < s cst H M= = e est + Ae est . (31) 

The way to use this variational principle is as follows. First, for a given J and U we 
compute the lowest eigenstate of H\ and this way obtain h. After the equivalence (14), 
finding the ground state of the local Hamiltonian Hi becomes equivalent to solving a 
modified Mathieu equation 

-U— -2J(n + l)cos(£)-2Jsin(£)— h = s est h, (32) 

which describes exactly the static properties of a pair of sites with open boundary 
conditions [28]. Once we have e es t we must still compute the correction Ae es t using a 
rather straightforward expansion which is shown in Sec. 3.4. Surprisingly, Ae est happens 
to be negative, so that it is actually an improvement over the simple estimate given by 
e est [See Fig. 5]. 

From the optimal variational state, — 0~ l \h)® M , and the estimate for the energy, 
£var = £est + A£ est , we may compute other observables. For the density fluctuations 
and nearest neighbor correlations we use the virial theorem 

t ^ 

( a j+l a j) = ~QJ £var ' 

<7 2 = ~^J £ var ~ n 2 , (34) 

whereas for other properties one has to evaluate numerically the matrix products shown 
in Sec. 3.4. This allows us to prove that for the product states |/i)® M correlations decay 
exponentially, opposite to what is expected in the supcrfluid phase, whose correlations 
should decay algebraically. Nevertheless, as we will see next, this family of states does 
estimate accurately the local properties of the optical lattice. 

3.3 Comparison to DMRG results 

We will now compare the results for the ground state energy, the correlation functions, 
and the variance of the particle number provided by the two variational ansatz, (18) 




U/J 

Fig. 3. (a) The ground state energy per site, e, (b) nearest neighbor correlation, 
ci = (cij +1 aj), and (c) variance of the number of atoms per site, a 2 = ((rtj — n) 2 ). 
Plots (b) and (c) use a log-log scale. The results of the DMRG (solid line) are 
obtained on a system with 128 sites, a maximum occupation number of 9 bosons 
per site and a reduced space of states of about 200 states. The estimates from 
the variational theory are plotted using dashed lines. The vertical lines mark the 
location of the phase transition according to [10]. The mean occupation numbers 
are denoted with circles (n = 1), diamonds (n = 2) and boxes (n = 3). 

and (32), and the Gutzwillcr ansatz [14] with those obtained by DMRG studies of the 
Bose-Hubbard model. The DMRG, developed 1992 by White [16, 17] in the area of 
condensed matter theory, is a very powerful numerical tool to investigate static and dy- 
namic properties of strongly correlated quasi-one-dimensional spin, fermionic or bosonic 
quantum systems. The DMRG is an essentially quasi-exact numerical method. The fun- 
damental ideas stem from real space renormalization methods: the system size is grown 
iteratively while the (exponentially diverging) size of the Hilbert space is kept constant 
by decimation. Hereby one tries to retain only that subset of states that is essential 
to describe the physical quantity under consideration. In DMRG these are expectation 
values with respect to low- lying states ( "target states" ) , and in particular with respect 
to the ground state wave function. 

DMRG builds up the system linearly: at each growth step, suitable density matrices 
for the target states are derived that yield information on the relevance of Hilbert space 
states. Building on this information, the states and operators are projected onto Hilbert 
subspaces of fixed dimension M containing the most relevant states. M is chosen to 
be small enough to be handled numerically, but large enough to obtain the desired 
accuracy; numerical results can be extrapolated in M to the exact limit of infinite M 
in the thermodynamic limit. However, results presented here have converged for the 
largest M considered and no further extrapolation was necessary. 



Details on the DMRG method may for example be found in [18]. In the case of 
the Bose-Hubbard model the DMRG has been used to study properties of the system 
[19, 20, 21]. The results of DMRG agree very well with exact diagonalization results for 
small systems, with quantum Monte-Carlo simulations e.g. [22, 23, 24, 25], and with 
13th order perturbation theory [26]. 

We have used the DMRG to study the properties of the ground state of the Bose- 
Hubbard model on one-dimensional lattices with 128 sites, and commensurate fillings 
n = 1,2, and 3. In Figs. 3(a-c) we show the results for the mean energy per site s, 
the nearest neighbor correlations, c\ = (aj +1 a :) ), and the variance a of the density, 
calculated both with the DMRG and with the variational estimates developed above. As 
expected, there are no indications of the phase transitions in these quantities, neither in 
the variational results nor in the numerical solutions. Rather, an inflexion of the nearest 
neighbor correlation points out the location of the superfluid-insulator transition which 
lies roughly between 3n and 4n (see [19] and ref. therein). The agreement of the two 
methods is fairly good above the phase transition and below it. 

A more detailed comparison is provided in Fig. 4 for the case n = 1. In this figure 
we plot together results from the DMRG, the variational ansatz derived above, the 
quantum rotor model and the well-known Gutzwiller ansatz. The Gutzwiller ansatz 
[27, 3] is a variational ansatz which reduces the wave function to a product of single-site 

wave functions, \^g) — ^IjLi where = X)m=o f™> \ m i) an< ^ /™ are constants. 
Such a wavefunction cannot be used in the one-dimensional Mott insulator regime, 
because a perturbative study of the Gutzwiller ansatz shows that the corrections of 
order 0(J/U) are lost and all correlations become zero. However this ansatz gives good 
results in the superfluid regime, where the long-range order is well described by I^g), 
and we can use these results and those of the DMRG to assert the accuracy of our 
variational estimates. As Fig. 4 shows, in the Mott insulator regime, the DMRG results 
agrees perfectly with our variational theory for the Bose-Hubbard model and for the 
quantum rotor model. Close to the phase transition is the point at which the quantum 
rotor model no longer describes well the atoms in the optical lattice due to the growth 
of density fluctuations. At this point we also observe a small disagreement between the 
DMRG and the coherent states, which is due to the growth of long range correlations 
and vanishes as we get deeper into the superfluid regime. 

3.4 Numerical evaluation of the upper bound 

In this section we will show how to compute the corrections to the local energy, Ae es t, 
from Eq. (31). We basically need a method to compute expectation values of the operator 
0~ 2 around product states which have the form 

\4>) = \h) 1 ■ ■ ■ \h) k ^\Ah) k \Bh) k+1 \h) k+2 ■ ■ ■ \h} M , (35) 

in which at most two contiguous vectors are affected by two single-connection operators. 
For instance, this is the case of H 2 \i/;), where A and B are E + , £~ or S z , and also of 
\ip), where A and B are just the identity. After some manipulations it is possible to 
write, for the optimal variational state \ij>) = -1 |ft)® M , 

M-2 

{0- 2 H 2 ) it = J2^ t °( HO ) k ~ lz '' ( HO ) M ~ k ~ 1 u, 
fc=i 

Z k = J(H z OH + -H_OH z )-U(H z OH z ), 
with the real matrices and vectors 
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Fig. 4. (a) The ground state energy per site, e, (b) nearest neighbor correlation, 
ci = (aj +1 aj), and (c) variance of the number of atoms per site, a 2 = ((nj — n) 2 }. 
Plot (b) and (c) are in log-log scale. Using filling factor n = 1, we show results 
from the variational model for the Bose-Hubbard model using phase coherent states 
(solid), the quantum rotor model (dashed), the Gutzwiller ansatz for the Bose- 
Hubbard Hamiltonian (dots) and DMRG (circles). Vertical dash-dot lines mark the 
location of the phase transition according to [10] . 
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We have used this technique to compute numerically the correction Ae es t using 
different lattice sizes and found small or no differences for more than 30 sites. Intuitively, 
this is because in the limit of large powers the matrices (HO) a become projectors on 
the eigenvector with the largest eigenvalue. 

Using the same type of expansion we may compute other correlations 



t O(HO) k - 1 {H-0) A {HO) 
tfO{HO) M u 



M-k+l . 



(36) 



For large values of A and large lattices, the numerator will decay exponentially as 7^, 
where 71 is the largest eigenvalue of the matrix H_0. 



Conclusions 



In this work we have studied analytically and numerically the properties of the ground 
state of an ensemble of bosonic atoms in an ID optical lattice. For the study of the 
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Fig. 5. The energy of the product ansatz contains a contribution from each con- 
nection, Eest, plus the interaction between neighbouring connections, Ae es t. In Fig. 
(a) we show that Ae e st (dash) is actually negative, and improves the estimate e ea t 
moving it towards the exact value, Edmrg (circles). Everything has been computed 
for n = 1. In Fig. (b) we show that the correction Ae es t does not change much for 
large lattices. 

atomic ensemble we have used both the quantum rotor model and the Bosc-Hubbard 
model. Exploiting the fact that in these models there exists only nearest neighbor hop- 
ping and local interactions, we have developed a variational wavcfunction that may be 
used to easily estimate local properties, such as the energy per well, the nearest neigbor 
correlations and the fluctuations of the density. In the case of the quantum rotor model 
we have verified our results with perturbativc calculations around the strongly inter- 
acting regime, and with a spin wave approximation around the superfluid regime. In 
the case of the Bosc-Hubbard model we have compared the variational estimates with 
numerical results obtained using the DMRG technique for a maximum density of three 
atoms per well. We have concluded that this procedure leads to fairly good estimates 
of local ground state properties of both Hamiltonians, in both the superfluid and the 
insulator regime, the largest disagreement being localized around the phase transition. 
In future work we forsee the application of these methods to time dependent problems. 
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